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We study vicinal crystal surfaces within the terrace-step-kink model on a discrete lattice. Includ- 
ing both a short-ranged attractive interaction and a long-ranged repulsive interaction arising from 
elastic forces, we discover a series of phases in which steps coalesce into bunches of nt steps each. 
The value of nj, varies with temperature and the ratio of short to long range interaction strengths. 
For bunches with large number of steps, we show that, at T = 0, our bunch phases correspond to 
the well known periodic groove structure first predicted by Marchenko. An extension to T > is 
developed. We propose that the bunch phases have been observed in very recent experiments on Si 
surfaces, and further experiments are suggested. Within the context of a mapping of the model to 
a system of bosons on a ID lattice, the bunch phases appear as quantum n-mers. 



I. INTRODUCTION 



The study of the equilibrium properties of a stepped vicinal crystal surface is important from technological and 
fundamental perspectives. The understanding of the surface morphology is key to phenomena like epitaxial growth, 
chemical etching and catalysis. In addition, it also provides a fascinating example of a problem which has a much 
broader context, reaching, via a mapping onto a one-dimensional quantum chain, into the realm of ID quantum 
liquids. In this paper we consider the morphology of vicinal surfaces in the case where step-step interactions are 
repulsive at large separations and attractive at short distances. Our work is aimed at explaining the features of the 
surface morphology of miscut Si (113) surfaces observed in two sets of experiments. The first is that that of the 
description of apparent tricritical phenomena observed in a beautiful set of experiments by Song et al. [0-^ on silicon 
surfaces oriented between the (113) and (114) crystalline directions. The second is the very recent observations of 
multiple-height steps by Sudoh et al. Q| on vicinal (113) surfaces misoriented towards a low symmetry azimuth which 
is directed 33° away from the [110] direction towards the [332] direction. 

Recent theoretical work |^J^ attempted to explain the experiments of Song et al. in terms of a continuum model of 
steps interacting via a long-ranged repulsive elastic interaction and a short-ranged attractive interaction. This model 
produces a tricritical point but requires somewhat artificial tuning to give the observed coexistence curve exponent and 
fails to describe the observed bunching of steps on the vicinal surfaces coexisting with the (113) facet. Furthermore, 
the experimental system seems to be very close to or in the region where the renormalization group method that is 
used to study the model fails. Here we explore the consequences of retaining the discrete, atomic nature of steps on 
a crystal surface within the same model. We discover entirely different physics. The steps do not phase separate but 
instead can coalesce into bunches whose size depends on the relative strengths of the short- and long-range interactions 
and temperature. Vicinal surface phases can be characterized by widely separated n-bunches and transitions occur 
between phases having different values of n. We propose that bunch phases with large n correspond to the two-phase 
coexistence region of Song et al and that the n=2,3,4 bunch phases produce the multiple-height steps seen by Sudoh 
et al. For bunches with large number of steps, we show that, at T = 0, our bunch phases correspond to the well 
known periodic groove structure predicted by Marchenko Q]. While Marchenko's calculations are valid only at T = 0, 
our theory predicts the evolution of the periodic groove structure as the temperature of the surface is raised. In 
addition, we show that Marchenko's phcnomenological parameters like the edge energy can be determined in terms 
of the parameters that characterize the long and short range step-step interaction strengths. The physics of bunching 
transitions that we study here is very different from the recently studied step bunching due to applied stress, electro- 
migration and other kinetic effects [ p^ . In our case, the bunching transitions occur in an equilibrium setting due 
to competing interactions acting on different length scales. The relevance of the bunching transitions described here 
go well beyond the crystal surface problem. In particular, the model provides opportunities to study many-body 
phenomenon arising from coupling of modes on different length and energy scales. A brief description of this work 
has already been reported ||] . 

The model we use is the terrace step kink model (TSK) [^. This lattice model can be mapped onto an equivalent 
quantum mechanical model of interacting spin-less fermions or hard-core bosons. Details are found in earlier work 
[p|-|ri|. In the quantum picture it is well-known that the 1-bunches form a Luttinger lattice liquid. Our results 
generalize this picture to a Luttinger lattice liquid which can form dimers (the 2-bunches), trimers (the 3-bunches), 
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and in general n-mers. The transitions between these phases are quantum many-body phenomena of a type hitherto 
unexplored. In the following section, we introduce the TSK model and develop the mapping to the hard core boson 
problem. In Section III, we analyze the model in what we call the extreme dilute limit (EDL), where bunches are 
widely separated. For T = 0, the exact solution clearly reveals the essence of our problem, a series of first-order 
bunching transitions. For bunches with large numbers of steps, we show the connection between the bunches phases 
and the periodic groove structure predicted by Marchenko 0. In the same EDL, using exact diagonalization on 
small systems to reveal the nature of the bunching transitions, we develop the physics for T > 0. Section IV uses 
phenomenology coupled with scaling arguments to include interactions between bunches in the case that they are 
widely separated and to study the effects of the interactions on transitions from one bunch size to another. In 
Section V we use Quantum Monte Carlo simulations to check the results of exact diagonalization by explicitly looking 
at particle-particle correlations through a series of bunching transitions. Section VI applies our results to crystal 
shapes. Section VII gives a detailed comparison of our results with continuum theories. In Sec. VIII we compare 
the key experimental findings with the predictions of the theory. Differences and similarities are highlighted, and 
directions for further development of the theory are suggested. We conclude by discussing the main results and 
pointing directions of future research. 



In order to study the thermodynamics of surfaces that are miscut by a small amount from a reference surface, we 
will first derive the surface free energy by treating the miscut angle and temperature as thermodynamic variables. The 
surface free energy will be expressed in terms of the ground state energy of a one-dimensional quantum mechanical 
model that will be derived using transfer matrix methods. In the study of interacting steps it is common to use both 
discrete and continuum models. In this section we introduce the TSK lattice model and the equivalent quantum 
mechanical model. We do not go into all the details of the derivation, but refer the reader to earlier work . We 
then compare the kink energy in the lattice model to the step-stiffness in the continuum model. It turns out that 
it is important to understand the limiting process which brings out the equivalence of the two models. While it is 
generally true that when the limiting process is carried out properly, it is possible to switch between the models, 
caution should be exercised against making erroneous conclusions about one model based on the results of the other. 
The bunching transition is one such example, where the cut off distance (lattice spacing) in the lattice model plays a 
crucial role. As we discuss in Sec. VII, there is nothing corresponding to bunching transitions in the continuum model; 
they are purely lattice effects. On the other hand, to study the interactions between steps whose average separation 
is much larger the lattice spacing, one can use either the lattice model or the continuum model. We will now turn our 
attention to the TSK model and (1) derive the quantum hamiltonian in terms of hard-core boson operators, (2) make 
the connection between the free energy of the surface to the ground state energy of the quantum mechanical system 
and (3) derive the connection between the kink energy in the lattice model and step stiffness in the continuum model. 



Let us consider a square lattice with M rows and L columns as shown in Fig. (p, where steps between terraces are 
specified by bonds on this lattice. Let the steps in the TSK model run parallel to the y-axis, so that the slope of the 
crystal face is in the a:-direction. One can specify a configuration of this system by specifying the a:-coordinates of 
each of the N steps at each of the y-coordinates. The lattice spacing in the x and y direction is denoted by ax and 
tty respectively. The density of steps is related to the slope of the interface under consideration through the relation 
s = tan{6) = {Nazi Lax) where is the lattice parameter in the z-direction or the step height. For notational 
convenience we specify a configuration by the state vector \xi,X2, ■..,XN)k, where Xj denotes the position of the jth 
step in the kth row. The position of the step can only change by ±1 between the rows by creating a kink. The transfer 
matrix describes the propagation of steps from row to row. It can be determined by a kink energy F, an energy rjo for 
a unit length of the straight step, and an interaction energy Vij between parallel straight pieces of steps labeled i and 
j. To render interactions between different rows (steps interacting at different values of y) tractable, we treat them 
in an average way: we compute Vij as the interaction between a unit length of step at i and an infinite straight step 
at j. We expect this approximation to preserve the qualitative physics of interest here. In going from row k to row 
fc -|- 1, if the ordinates of p steps are unchanged while the other (iV — p) steps change by ±1, then the transfer matrix 
element has a contribution exp {—fSayNrjo — (3T{N — p)), where /3 = l/ksT. The matrix element vanishes if Xi — Xj 
for i 7^ J, i.e. the steps do not cross. The transfer matrix element also has a multiplicative element contribution from 

the step interactions, that can be written as exp ( —Pay X]i<7 ) • ^ convenient operator representation is obtained 



II. HAMILTONIAN OF INTERACTING STEPS 



A. Terrace Step Kink Model 
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by associating a step creation operator Oj for a step at xj in row k: OjCij denotes a straight step while aj+iflj and 

a]aj+i denote a kink to the right and to the left, respectively. The a's are defined so that (a])^ = 0, which imposes 
the non-crossing condition. This gives us a choice of a fermionic or hard-core bosonic representation of a's. We adopt 
the hard-bosonic representation, where a's satisfy the usual bosonic commutation relations 



4, a] 



ai, a 



= 



when i ^ j and satisfy anti-commutation relations on the same site, i.e., 



tti, a. 



4>4 



0, 



1. 



(1) 



(2) 



The onsite anti-commutation relations take care of the hardcore condition (a^)^ = 0. This representation is equivalent 
to the fermion representation introduced in Ref. |^ . The fermion operators that anti-commute on different sites are 
related to the hard core boson operators through the unitary Jordan- Wigner transformation ||l3| . 

Let us denote the (N x N) transfer matrix by T. For mathematical convenience we take the continuum limit in the 
y-direction (i.e. we take the limit of the lattice spacing Uy —^ 0). In this limit to first order in Uy the transfer matrix 
can be written as 



r = 1 - 



knT 



where 1 is the identity operator and the quantum hamiltonian is defined as 

H = r]{T) ^ - ^ ^ [al+iOi + ai+ia| - 2n.jj + ^ Vij 



(3) 



(4) 



where {i,j) is used to denote a sum over distinct pairs of bosons, rij — aja^, and the free energy of a straight step 
TjiT) — rjQ — 2t, where the hopping parameter t is defined as 



lim 



exp i-PT) 
f3ay 



(5) 



It should be noted that the kink energy F is a function of ay and diverges in the limit of vanishing ay, but the ratio 
given by Eq. (|^) remains finite. Note that since t is a monotonically increasing function of T which vanishes as T — > 0, 
t itself serves as an effective temperature variable. 

The eigenvalues of the transfer matrix can be related to the projected surface energy of the miscut surface through 
the relation [B| 



lis) 

cos{e) 



70 



keT 
MLaxay 



log (Tr f^' ) 



(6) 



where 70 is the surface tension of the reference surface and lis) is the free energy per unit area of the miscut surface 
i.e. the surface tension. In the continuum limit where Eq. (|^) is satisfied, the projected free energy can be related to 
the ground state energy Eq [N, T] of H through 



/(s) = 70 + 



Ea[N,T] 
Lax 



(7) 



This relation can be used to relate the physical properties of the surface like the curvature, stiffness etc. to the ground 
state energy of the hard-core boson system. We now turn to a discussion of the form of the step-step interactions 
that will enable us discuss the physics of faceting transitions. 



B. Step-Step Interactions 



In this paper we consider step-step interactions consisting of a repulsive long range part and an attractive short 
range part. We write the inverse square long range part in the form 
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and the short range part as 

KI = -E^'=V.l.fe ' (9) 
fe=l 

where 5 is the usual kronecker delta. The quantity G has contributions arising from the elastic interactions as well 
as from the electronic charge rearrangements along the steps . The elastic part depends on the elastic constants of 
the crystal; for isotropic crystals it has the form p4| , p^ 

where E and a are the Young's modulus and Poisson's ratio respectively and / is the strength of the vector surface 
dipole force at a step. In general for non-metallic crystals the electronic rearrangements along the steps can be 
ignored, and G > 0. In the case of metallic crystals, G can have either sign depending on the details of the electronic 
charge rearrangement along the steps. The short range potential depends on the details of the step-step interactions 
on atomic scales and has to be worked out from a microscopic calculation. There is ample experimental evidence 
[p|-p|, p^|jT^ ] that this quantity can be attractive, i.e. Uk > 0, on a few lattice spacings from the step iQ. The purpose 
of this paper is to explore the phases of vicinal surfaces and the associated crystal shapes in crystals with G > 
and Uk > 0. In the remainder of the paper we will consider the case = 1 and Ui = U, though in Sec. VII we 
comment on the physics of systems with ris > 1. For convenience, we re-scale the hopping strength to unity, so that 
the interaction parameters get scaled to G ^ g = G/t and U ^ u = U /t, while the free energy per step is rescaled as 
f]{T) — •q{T)/t. From this point on, we will focus on the rescaled Hamiltonian TL, written as 



1 (u 



rurij. (11) 



Note that g scales inversely with the effective temperature variable t. 

C. Step Stiffness and Kink Energy 

In the study of the thermodynamics of steps it is very common to use a continuum model, where the Hamiltonian 
of a single step with coordinate y{x) is written as 

"-^[('1)'^.. (12) 



2 J \dx, 

where E is the step stiffness. For a single step it is instructive to make the connection between the kink energy F 
in the lattice model and the stiffness E in the continuum model. To this end, in Appendix A we derive, within the 
lattice model, the free energy of the step. By comparing the energies of the continuum model and the lattice model, 
we see that the the step stiffness can be written as 

E^hm^f (13) 
a„^o 2a2 Vexp(-/3F)y 2alt ^ ' 

It should be noted that the quantity in the brackets assumes a finite value in in the limit of vanishing ay . The above 
equation is key in understanding the connection between the discrete and continuum models. We will return to this 
equation when we discuss the stiffness of step bunches. 

III. PHASE DIAGRAM IN THE EXTREME DILUTE LIMIT 

The form of the interaction that we have considered for the step interactions has competing components acting 
at different length scales. As noted before, we will focus on the case where there is an attractive potential on the 
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near neighbor site only i.e ris = 1 and Ui = U. We now show that a number of features of this many body system 
can be inferred from exact diagonahzation of small systems in conjunction with some simple analytical calculations. 
In the discussions to follow, we will present simple physical arguments to show that the for sufficiently attractive 
interactions, the steps on the surface rearrange themselves into bunches at low temperatures. At T = 0, where the 
entropic effects are unimportant, the bunch size is completely determined by the ratio of the strengths of the attractive 
and repulsive interactions i.e. U /G. With increasing temperature, the steps start peeling off from the bunches in 
a series of bunching transitions. In the extreme dilute limit (the details of which will be clear in the discussions to 
follow), we obtain a phase diagram that shows these bunching transitions as function of temperature. 



A. Phase diagram at T = 



Let us now consider the limit u — > oo and g — > oo but u/g finite, so that i — *■ 0, and the hopping part of the 
Hamiltonian Eq. ( pj] ) can be completely ignored. This corresponds to taking the zero temperature limit of the vicinal 
surface. In this limit, the energy of the system is obtained by minimizing the total potential energy, and depends on 
the single parameter U/G. A little reflection reveals that in the dilute limit, the minimum energy configuration for a 
given U/G consists of bosons in bunches of size nb that are well separated from each other. In a bunch, the bosons 
sit next to one another. For a bunch of size ni,, the energy per boson is given by 



e('^&) = = '7(0) - U 



(«6 - 1) 



G 



«(,— 1 



1 (rtf, - i) 



OiiN/Lf), 



(14) 



where the term 0{{N/L)'^) comes from bunch-bunch interactions and is small in the dilute limit. We call the limit in 
which the bunch interactions are completely ignored, the "extreme dilute limit "(EiDh). A straightforward minimization 
of e{nb) in the EDL reveals the phase diagram shown in Fig. (|^), where the circles marked on the U/G axis separate 
bunches that differ in size by one. With increasing U/G, the bunch size keeps increasing; for U/G < I the bunch size 
is 1, for 1 < U/G < 1.5 the bunch size is 2, for 1.5 < U/G < 1.83 the bunch size is 3, for 1.83 < U/G < 2.08 the 
bunch size is 4 and so on. At the circles marked in Fig. (||) (i.e. at U/G = 1, 1.5, 1.83, 2.08, ...) there is a coexistence 
of bunches that differ in size by one. From the graphs of the energy per boson vs. U/G shown in Fig. (||) for different 
rib's, we see that the slope at the points where these curves intersect are not the same. The T = bunching transitions 
brought about by changing U/G can then be considered first order. 

With increasing U/G, the spacing between the regions that differ in bunch size by one becomes smaller and smaller. 
For large U/G, Eq. (14) becomes 



e G 1 

e{nb) = eoo + — Inn^ + O(^), 

rib rib f^f 



(15) 



where too = ^7(0) + GC(2) — C/ is the energy per step in a large bunch, = U — G — GC is the energy associated with 
the ends of the bunch, and C = 0.577 is Euler's constant. Minimizing e{nb) in Eq. ( p^ ) with respect to rib leads to a 
bunch size given by 



rib « 0.561 exp(C//G), C//G > 1, 



(16) 



which indicates that bunch size diverges rapidly with increasing U /G. 

The above simple calculation has shown that the steps split themselves into bunches in the EDL at T = 0. In the 
next subsection, we study the unbinding of the steps with increasing temperature in the EDL. First, however, we 
connect the result Eq. (16) for large bunches with earlier work of Marchenko |Q. 

For the case at hand, if elastic interactions are ignored, steps have only the short-ranged attractive interaction ~U. 
Here there is a coexistence (see, e.g. |^) between the flat, step-free facet and the facet formed when there are steps at 
each site. The latter facet is inclined at an angle 6*0 = tan~^{az/ax) with respect to the former. All facets with angles 
9 such that < 9 < 9q are unstable and will decompose into coexisting facets with angles and 6*0. Marchenko 0| 
showed that when elastic interactions are taken into account, facets with < 9 < 9o become stable, being formed from 
an array of grooves as shown in Fig. (^). Each groove is composed of a step-free portion followed by a step bunch, 
in the terminology of this paper. In fact this "groove" phase looks just like a bunch phase with rib large. It should 
then be possible to connect Marchenko's results with those of Eq. (|l5| ) in the limit, the EDL, where the bunches are 
widely separated. This corresponds to the case of small 9. Marchenko also assumes small 6*0, as will we in making 
the connection to his work. In terms of the geometry shown in Fig. (M), we have 
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La sin0 LgO , , 



where Lg is the length of a Marchenko groove and a ~ ^ 0^^ + a^. The surface energy per unit area of the grooved 
surface is 

^ ^ 70^2 + 7(^0)-^! 

^9 

_ 70(^2 +il COS 6*0(1 + (£00/7003;)) G Yj^(^^9^ \ 



Lg Lg 0*6*0' 

Lg "''a*9o' 



70(e) -^ln(^), (18) 



where 



and 



7(eo)=7ocos^?o + ^^, (19) 



a* = a^e"^/^. (20) 
Eq. ( |l8| ) is precisely Marchenko's result in the EDL if G can be identified as 

in which F is the surface force acting at the ends of the bunch. Now, an expression for G was given in Eq. (|l^) in 
terms of /, the strength of the vector surface dipole force at a step. To connect F with /, note that the surface force 
density due to an array of steps with density n(x) is given by 

f{x) = J fS'ix - xoHxo)dxo = /^^. (22) 

J-{x) can be integrated to obtain the force F on the surface at the edge of a bunch. Noting that the step density has 
a discontinuity of size no, where no — (l/ax) is the uniform step density in a bunch, we obtain 

F^ff'-^dx^l. (23) 
J dx 

This result proves the equivalence of Eq. ( pT| ) to Eq. ([l^) , and demonstrates that our theory is in fact a microscopic 
version of Marchenko's. In addition to connecting the surface force to the strength of the force dipole at the step, 
Eq. (^o|) gives an microscopic interpretation to the phenomenological edge energy, eg, in the Marchenko theory. 
Explicitly, in terms of the microscopic step interaction strengths, ee = C/ — 0.4236". 
Minimizing 7^ with respect to Lg at fixed 9 and ^o gives 

L9 = (24) 
which is Marchenko's result in the EDL and is equivalent to Eq. (|6|). At this minimum, we have 

lM^lo[9)--^. (25) 
ea*9Q 

To further understand the grooved configuration, we make use of the relations Lg sin 9 = Li sin 9q and Lg sin(0o — 
9) = Li sin 6*0 to show that 

7o(e)=cos(0)[7o + ^^^]. (26) 
tan(6^o) 
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The full projected free energy is, for small 6 and Oq, given by 



/.W = ^=7o + (^-^)f. (27) 
cos 61 fla; ea* 6'o 



The contribution to the projected free energy at 0{9^) comes from bunch-bunch interactions that are ignored in the 
EDL. In the limit of vanishing 9 that we are considering, the bunches interact via an inverse square interaction of 
strength Gn^. A straightforward computation of the bunch interaction contribution yields 

igi.9) _ , ,eoo G .e Gtt^ (e^ 



cosfc* ax ea* Oq 6ea* V^g, 



Marchenko's full result replaces Eq. (25) by 



7s(e)=7o(e)-^sin(7r^). (29) 
7ra uq 



When expressed as an expansion in 0, up to 0{6^), for small ^o, Eq. ( P9| ) becomes identical to Eq. (p8|), once again 
confirming that our model is identical to Marchenko's model at T = 0. In the following section we derive the finite 
temperature extension of Eq. ( |2^ ) . 

B. Phase Diagram for T > 

At finite temperatures, it becomes favorable to form kinks in the steps for entropic reasons. This causes the steps 
within a bunch to fluctuate from their mean positions, leading to a decrease in the free energy. On general physical 
grounds then, we expect the free energy per step in larger bunch sizes to decrease less rapidly compared to smaller 
bunches as steps within larger bunches will have their fluctuations more constrained. If this happens, the free energies 
of smaller bunches will fall below those of larger bunches. We then expect larger bunches to rearrange themselves 
into smaller bunches, with increasing temperature, until eventually at very large temperatures only one step is left in 
the bunch. In this section, we show that this indeed happens. 

In order to compute the bunch size at a given temperature in the EDL, it is not necessary to solve the step problem 
in the thermodynamic limit, or equivalently the hard core boson for a large number of bosons confined to very 
large number of sites. Instead, computing the energy per boson in a bunch by solving the bunch problem by exact 
diagonalization is sufficient. This simplification can be understood by noting that bunch-bunch interactions give rise 
to corrections to the ground state energy that are smaller than the leading bunch energy by a factor s^, where s is 
the boson density introduced earlier. As a result, in the limit that the interaction between the bunches are ignored, 
comparing the energy per boson in bunches of different sizes will determine the stable bunch size. Even computing 
the ground state energy of a bunch confined on a large number of lattice sites could be difficult due to demands on 
computational resources. However, for studying the transitions from bunch of size rif, to rifc — 1 (rif, > 2), the bunches 
are sufficiently bound that a small number of lattice sites (typically twice to thrice the bunch size) is sufficient to get 
accurate answers. We will present more details on this in the discussions to follow. 

Following the notation introduced for T = 0, we write the energy per boson as 

E^\M2l^ = ^(D + fin,, T) + 0{{N/Ln (30) 

where in the EDL, the bunch size U}, is obtained by minimizing f{nh,T), which is the contribution to the energy 
arising from interactions within a bunch. It is obtained by solving for the ground state energy of Hamiltonian given 
in Eq. (|ll|) for ni, bosons. We illustrate the results of this procedure in Fig. (||) for u/g = 1.65. Here the energy per 
boson is plotted as a function of inverse temperature g = G/t for ni, — 2,3 and 4. At low temperatures {g > 6.8), 
the energy per boson in the 3-bunch the lowest value, while at moderate temperatures (6.8 > g > 4.54) the 2-bunch 
gives lowest energy per boson, while at high temperatures {g < 4.54) the 1-bunch has the lowest energy (the 1-bunch 
becomes stable when /(2,T) = 0). From our results we have computed phase boundaries, shown as the solid curves 
in Fig. (^), for transitions from 1-bunch to 2-bunch, from 2-bunch to 3-bunch, and from 3-bunch to 4-bunch. In order 
to obtain the phase boundaries between the n^-bunch phase and the ni, — 1 phase for a fixed u/g, one has to find the 
value of g that satisfies the equation /(rib, T) — f{ni, — 1,T) — 0. We do this numerically by using a standard root 
finding algorithm which approaches the root iteratively ]l9| . The root finding algorithm requires repeated evaluation 
of /(rib, T) — f{ni, — 1, T) using exact diagonalization of the nt and ni, — 1 bunches for different values of g. 
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A key feature of plots shown in Fig. (||) is that the energy curves for bunches that differ in size by one intersect with 
a finite slope. This means that the derivative of the free energy of the vicinal surface with respect to temperature, 
or the entropy, is discontinuous at the transition, implying that the the transitions are first order in the EDL. In the 
next section we will argue that, when bunch- bunch interactions are taken into account, the discontinuity in slope can 
vanish and the transition becomes continuous. 

In calculating the energy of the 3-bunch and 4-bunch, we confined the bosons to 12 sites on a ring, while for the 
2-bunch, 200 sites. The latter is equivalent to a 1-body problem describing the relative motion. These ring sizes are 
sufficient to get accurate answers because the bunches are tightly bound with a mean particle spacing between 1 and 
1.5, so that as long as the ring is bigger than about twice the bunch size, the energies do not change much. We have 
also computed the ground state energies using the GFMC method for bunches confined on lattices with up to 100 
sites and verified that the energies plotted in Fig. (||) are accurate. In Fig. we show the energies computed for 
a 3-bunch on 12 sites using exact diagonalization and on 30 and 100 sites using GFMC( see Sec. VI for details on 
GFMC simulations). For sake of completeness, we have computed the mean particle spacing for the 3-bunch along the 
phase boundary separating the 2-bunch and 3-bunch regions and for the 4-bunch along the phase boundary marked 
between the 3-bunch and 4-bunch regions (see Fig. (^)). These spacings vary from about 1.3 to 1.6 lattice spacings 
as u/g is varied from its values at the transition onsets up to u/g — 3. It is clear that the bunches are very tightly 
bound in all the cases. 

It is of considerable interest to understand how the limit 5 — s- is approached for our model. This limit, where 
only a near-neighbor interaction remains, has an exact solution (see, e.g. Refs. ]9|,|20|]) involving a first-order phase 
transition from a gas (infinitely dilute) of 1-bunches to a condensed phase (an 00-bunch). Note from Fig. (^) that 
the 2-3 boundary crosses the 1-2 boundary at u/g = 2.4 so that a 2-bunch is no longer stable for u/g > 2.4. We also 
show the 1-3 transition line extending from u/g = 2 upwards. This line cannot be distinguished from the 2-3 line 
for u/g > 2.4, but in fact lies slightly below it and above the 1-2 line. Thus, for u/g > 2.4 in Fig. (||) the 1-bunch 
goes directly to the 3-bunch and then subsequently to the 4-bunch as the temperature or 1/g is lowered at fixed u/g. 
Further computation for values of u/g up to 8, depicted in Fig. (||), reveals that the 3-4 line comes very close to the 
1-3 line while remaining below it. We then conjecture that as g — > the n-1 transition lines become asymptotically 
equal for arbitrarily large n. This is clearly an approach to the l-oo transition, which is the first-order transition in 
the pure short-ranged potential case. As the transition in the short-ranged case is at u = 2 the asymptotic 

slope of the transition lines should be at {\/g)/{u/g) = 1/2. This is indeed the case as seen in Fig. (||). 

IV. PHASE DIAGRAM IN THE DILUTE LIMIT 

We will now consider the effect of the bunch-bunch interactions that were ignored in the previous section. We 
adopt the following steps in discussing this effect: (1) Using the general principles of quantum mechanics, we argue 
that when the the inter-bunch interactions are turned on, the first order transition should become continuous. (2) We 
then estimate the regions around the first order transition lines (obtained in the previous section) where the inter- 
bunch interactions become important. (3) Finally, we develop a perturbative approach to compute the free energy 
to 0{N/L)^ in the region where the inter-bunch interactions are small. In the next section, we will then present the 
results of GFMC simulations on larger systems and show them to be consistent with our simple theory. In particular 
we compute the 2-point correlation function, which clearly shows the bunching transitions in large systems. 

The effects of interactions on the bunches can be understood qualitatively using arguments based on general 
principles of quantum mechanics. We confine our discussion to those regions in the phase diagram where the nt-nt — 1 
phase boundary is well separated from all the other phase boundaries. Consider Fig. (|^), where we schematically plot 
the energy per step in bunches of size nt and rib — 1 as a function of the inverse effective temperature g for some fixed 
value of u/g. This is similar to Fig. (|^), which was used to study phase transitions in the extreme dilute limit. In this 
figure, we identify three regimes, namely, (1) regions where the energy per boson in the nf,-bunch is much lower than 
the Uf, — 1-bunch, (2) regions where the energy per boson in both the Ub-bunch and the Ub — 1-bunch have comparable 
values and (3) the regions where the energy of the nb — 1-bunch is much lower than the nb bunch. In regions (1) and 
(3) the ground state energy of a system with large number of bosons can be estimated using perturbative methods. 
Here to leading order the energy of the system is the energy of a bunch times the total number of bunches. The 
next to leading order comes from the bunch-bunch interactions via an inverse square interaction. In region (2), called 
the "critical region" , the physics is more complicated, since the system has the nearly same energy in both the Ub 
and rif, — 1 bunch phases. In this case we expect a splitting of the energy levels as sketched in Fig. (||), where the 
degeneracy of the energy levels is shown to be lifted due to bunch-bunch interactions. The extent of the region with 
significant splitting is denoted by Ag{nb ni, — 1). This width, which vanishes in the limit of vanishing densities, 
can be written as 
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Ag{nb ^rib-l)-^ AT,{nb -> rib - 1) ~ 5^"^ , (31) 

where ATdrib Uh — l) is the width, in temperature, of the critical region and the exponent > wiU be computed 
later in this section by comparing the free energies of the bunches of size rif, and rib — 1, derived perturbatively. The 
physics in the critical region is of general interest, and its relevance goes well beyond the problem we are studying. It 
involves coupling of modes with very different length scales. A tightly bound n;,-bunch in region (3) can be thought of 
as having phonon excitations involving all the bunches as well as excitations within the bunch. While the wavelength 
of the former type of excitations is large (of the order of bunch spacing s~^), the latter modes are localized in regions 
of the order of few lattice spacings. A bunch of size Ub — 1 in region (1) will also possess similar excitations on both 
large and small scales. While the couplings between the excitations on different scales are small in regions (1) and 
(3), we expect strong coupling between the short and long wavelength modes in the critical region. 



A. Perturbative Treatment for Stable Bunches 



In order to compute the bunch-bunch interaction energy in regions with stable bunch sizes (like regions (1) and (3) 
in Fig. (^), we make the following assumptions: (1) If a bunch has Ub steps, then the stiffness of a bunch is Sb = ribS, 
where S is the single bunch stiffness in the continuum limit. (In the discrete model, the kink energy required to make 
a kink in the bunch becomes UbT. However, taking the appropriate limit of vanishing ay would guarantee that the 
combination exjp{— (3nbT) / ay would scale like n^). In making this approximation, we have ignored the contribution to 
the bunch stiffness that arises from interactions of the steps in the bunch. This is justified by the fact that average 
spacing in a bunch, for u/g < 3, remains well under two lattice spacings, so that the steps in a bunch remain tightly 
bound as a single entity. The approximation will be less accurate if the step spacing in the bunch is large, in which 
case bending of steps would also alter the interaction energy within the bunch. (2) The bunches interact with each 
other with a renormalized inverse square potential of strength = nig. This approximation is justified by noting 
that, in the dilute limit, the spacing between the bunches is very large, so that one bunch sees other bunches as if 
they were single lines without any internal structure. We note that the contributions to the energy of the system 
arising from the short range bunch-bunch interactions are smaller than the inverse square contribution by a factor 
proportional to the density of steps s ^ 1. This can be shown very easily using a Hartree-Fock estimation of the 
short-ranged contribution pl[ |. The energy of the system can then be computed using the Calagero-Sutherland model 
of hard-core bosons interacting via an inverse square law, for which an exact solution was provided by Sutherland 
p2[ . Using his solution, we can write the energy per site in a system with bunches of size Ub as 



Eo[s,T,nb] [r]{T) + f{nb,T)]a^ 



Gub'^afYi 



0(3% 



(32) 



where 



Ab = ^ 



l+Jl + 2nlg 



(33) 



Using the above form of energy, one can compute such quantities as surface stiffness, crystal shapes etc, which we 
will turn to following a discussion of the physics in the critical region. An important point to make at this juncture 
is that in the limit of zero temperature, for large Ub, Eq. (|3|) is identical to Marchenko's result given in Eq. ([2^). It 
is then clear that Eq. (^) is the finite temperature extension of Marchenko's result. 



B. Analysis of the Critical Region 

In the critical region, since the energies of both the Ub and rib — 1 bunches are nearly equal, we retain the terms 
of order s'^ arising from bunch-bunch interactions. In order to estimate the width of the critical region, we use the 
criterion E{s, T, nt, — 1) w E{s, T, rib), to obtain 

s^^\f{nb,T)-f{nb-l,T)\. (34) 

If 

|/(nb,T) - /(71b - l,T)\ ^ \Ag{nb rib - 1)1""', (35) 
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we obtain the result =2/ , for introduced in Eq. ( |3l| ) . For all the transitions from the 4 to 3 bunches as 
well as for the transitions from the 3 to 2 bunches we see that the curves for the energy per step intersect with a finite 
slope for the parameter range that we have studied. This means that for these first order transitions — — 1, 
the width of the critical region scales like 

ATc (X (36) 

For the unbinding of the 2-bunch to the 1-bunch, we find, for 5 < 3/2, 02 = 2/-yi + 2g, while for g > 3/2, a2 — 1, 
from which one can infer that 

ATc{2 ^ 1) oc s^'''^, g < 3/2 

cx s^ .g>3/2. (37) 

As noted earlier the 3 bunch to 1 bunch transition preempts the 2-bunch to 1-bunch transition for g < 1.86. The 3 
bunch to 1 bunch transition was seen to be first order in the entire range investigated. 

Note next that the simplest conjecture for the shift of the transition temperature due to interactions is that it is of 
the same order of magnitude as the level shift given in Eq. (|l]). Explicitly, for the transitions that are first order in 
the EDL, 

\Tc{nb n[; s) - Tc{nh n[; 0)| cx . (38) 

While we find that all the transition in the model are first order in the EDL, it is interesting to compute the shift in 
transition temperature for the 2 bunch to 1-bunch transition temperature when g < 3/2. With our conjecture, and 
using Eq. ( p7| ) for g < 3/2, one obtains 

|T,(2^ l;s)-r,(2^ 1;0)| ocs^^. (39) 

Note that for this case our conjecture indeed yields the result derived within continuum theory |^,^ using renormal- 
ization group methods. 

While we have used simple arguments to predict the width of the critical region and the shift of the transition 
temperature, a detailed analysis of the nature of unbinding transitions and the form of the free energy require further 
study of the many body problem involving the thermodynamic limit. One natural way to study this is by the means 
of computational methods like GFMC or exact diagonalization methods for large systems. In the following section, 
we discuss the results obtained using GFMC simulations. 



V. QUANTUM MONTE CARLO SIMULATIONS 

In the preceding section we argued, on the basis of results for the extreme dilute limit from exact diagonalization 
of small systems, that the bosons undergo bunching transitions as strengths of the interaction parameters are varied. 
In this section, we will perform calculations on larger systems to directly study the dilute limit and attempt to 
understand the physics in the "thermodynamic limit" . In addition to exact ground state energy, we also study pair 
correlation functions in order to probe the structure of the many-body state, in particular the occurrence of bunching 
and transitions between different bunch sizes. 

The primary method we use is Green's function Monte Carlo (GFMC) a Monte Carlo method to study 

ground state properties of quantum many-body systems which is exact for bosons. Variational Monte Carlo (VMC) 
calculations are performed prior to GFMC, in order to optimize our variational wave function, which is used to 
generate the initial state and also as the importance function in GFMC. In Appendix II, we present a brief self 
contained treatment of these methods. The algorithm we use follows closely that of Refs [ p4 
reader to for more details. 

The variational many-body wave function used is of the Jastrow form, written as a product of two body correlations, 
i.e., 

where I the mean distance between particle i and particle j. The function /(r), where r is a positive integer, 

is chosen to be of the form 



25 1, which we refer the 
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f{r) = a{r - r^Y + r < ro 

= + r>ro, (41) 



Defining /(I) = /i, we choose = {a + ro)r^-\l-d), b = l + a/ro, b = af2e''^^' /ibr^+^) and a = {-d + h)/{l~ro)^ 
thus reducing our variational parameters to four, namely, /i, tq, d and a. We note that the parameters are chosen 
such that f(r) and df / dr at r = rg are continuous. These parameters have direct physical implication as shown in 
Fig. (0). 

In both the VMC and GFMC, we used 500 walkers (on the average for GFMC) for computing the averages. 
In the variational calculation, measurements were made after 10,000 metropolis steps. In the GFMC calculations 
measurements are made after 10,0000 steps during which system relaxes from the initial variational distribution. 
Measurements were made for a total of 20,0000 steps, and statistical errors were obtained by dividing these into 200 
blocks of 1000 steps each and estimating the variance of the block measurements. 

We present results for three separate cases, namely, u/g = 1.25, 1.65,2.0. These cases are chosen so that in the 
limit of 5 — !■ cx) the ground states have bunches of size 2, 3 and 4 bosons, respectively. 

First we show in Table I a comparison of GFMC and exact diagonalization results for the energy of a small system 
of 3 bosons. In the exact diagonalization the bosons are confined on 12 sites in a ring, i.e., with periodic boundary 
condition. The results are shown together with those from GFMC for the same system. The agreement is exact. This 
is to be expected, since GFMC yields exact ground state energies. Also shown are the energies computed by GFMC 
for a ring of size 30 and 100 respectively. We see that they are indistinguishable from the corresponding results for 12 
sites, indicating that at 12 sites the finite-size effect is already negligible. This confirms our earlier conclusion that the 
3 bosons are tightly bound. In Fig. (Q), we again compare GFMC results for 2- and 3-boson systems on 30 sites with 
our exact diagonalization results, over a wider range of parameter space. The excellent agreement reassures that our 
GFMC code is well-behaved, and that the small system sizes used in exact diagonalization calculations of the EDL 
are justified. 

In table II, III, and IV we present the variational wave- function and the energies Evmc and Eqpmc and compare 
them with the approximate energy Eapprox of the rescaled Hamiltonian given by Eg. (pi]) . (For sake of clarity we do 
not include the first term which is proportional to the number of bosons, since it does not directly affect the physics 
of bunching transitions.) For N bosons on L sites, Eapprox is computed from the energy per boson f(ni,,T) computed 
using exact diagonalization and the approximate interaction energy using 



^approx - + 3„4^3 ' ^^^^ 



where Af, is given by Eg. (p3|) . In computing Eapprox, the value of rib is obtained from the the phase diagram given 
in Fig. (^. For all the values of u/g, for large g, we see that there is good agreement between the GFMC and 
approximate energies and fair agreement for small values of g. These results confirm the approximations made in the 
analysis of bunching and bunch-bunch interactions. As g is decreased, the system is more guantum mechanical and 
each bunch is less tightly bound, hence Eapprox becomes worse. Note that the in all cases Eqpmc < Eapprox, which 
is not surprising since we did not consider the role of the short-range attraction between the bunches in evaluating 



Eapprox • 



We also calculate the pair correlation function which for L even is defined as 

L 
2iV2 



5(^) = :j^<II'^(2^-2;, +x,0>, (43) 



for l<a;<_L/2 — 1. The distance is halved due to periodic boundary condition. In the expression for g{x), Xi and 
Xii are the coordinates of the bosons i and i' . 

In GFMC the mixed estimate is used to compute the average indicated in Eg. (^3|) . Recall that the mixed estimate 
[ p3[ for the expectation value of an operator O is 

_ (^tIOI^o > .... 

- <Vt|Vo>' ^ ^ 

where is the ground state wave function while |^t) is the trial wave function. While exact for the energy, 
this estimator yields an approximate ground state expectation value for other operators. If \4't) is good, a simple 
extrapolation |23| ] of the mixed estimate with the variational expectation can lead to a further improved result. Our 
IV't) is relatively poor, giving variational estimates of g{x) which are rather structure-less in long range, so we did not 
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attempt this extrapolation. However, as we see below, even with this |'0t), 5(2;) 's computed from the mixed estimate 
clearly shows the bunching structure. 

For u/g — 1.25, Fig. (^l]) shows the pair correlation function for 10 bosons on 100 sites respectively. We see that for 
large g (g = 20 and g = 15), the pair correlation function shows two equal peaks near x — 20 and x = 40. Taking into 
account periodic boundary condition, we see that this indicates the presence of 5 2-bunches in the system. Note that 
the 2-bunch at 5 = 15 is more loosely bound than the 2-bunch when g = 20. When g = 10 the system is a completely 
unbound state. The case 17 = 12 is in the intermediate region, where the pair correlation still has a large peak at the 
nearest neighbor site, but does not show peaks at positions which indicate the presence of a 2-bunch. This value of 
g corresponds to a value of the potential such that the system is in the critical region of the phase diagram. In this 
region the 2-bunch ground state continuously unbinds into one with single bunches. 

We next consider the case u/g = 1.65, where we study the ground state of 12 bosons on 120 sites. Fig. ( p^ shows 
that for large g {g ~ 15 and g = 7.5) there is clearly a 3-bunch, which starts unbinding at 17 = 6. At g = 5 there is 
a 2-bunch. This starts to unbind around g = 4 to a single bunch ground state. It is clear that the cases g = Q and 
g — 4: fall in the critical region, where the unbinding of a 3 to 2-bunch ground state and 2 to 1 bunch ground state 
takes place, respectively. From Fig. (^3|), for g ~ 2, it is evident that a 4 bunch first unbinds to a 3 bunch, then to a 
2-bunch and finally a 1-bunch. 

In Fig. (|l^, we show points corresponding to Figs. ([Tl|)-([r^) on the phase diagram derived for the extreme dilute 
limit (EDL). Since these points are for larger systems in the dilute limit, with multiple bunches present and interacting, 
the precise phase boundary locations for these systems are expected to differ from those for the EDL. We see that 
the behavior of the pairing correlation functions is completely consistent with the conjectured phase diagrams. 



VI. EQUILIBRIUM CRYSTAL SHAPES 

In this section we will discuss the equilibrium shapes of crystals that have step-step interactions with competing 
components. The crystal shapes consists of both smooth and rough regions. The latter regions consists of step bunches 
separated from each other by distances depending on the local curvature of the surface. From the discussions on the 
phases of vicinal surfaces, it is clear that the size of the bunches depend on the temperature of the crystal. If the 
bunches have size, say rif,, at low temperatures, the size progressively decreases as the temperature is increased. In 
the temperature regimes where we can apply perturbative treatments to compute the step energies, we will find the 
crystal shapes for both cylindrical and 3-dimcnsional crystals. 



A. Cylindrical Crystals 

We start with the simple case of a cylindrical crystal where the variations in the shape occur in a 2-dimensional 
plane i.e, the crystal is infinitely long in the third direction. The crystal shape is obtained by using the projected 
free energy of a surface of orientation 9 with respect to the reference surface using Wulff construction. If one tunes 
the temperatures such that bunches of size rif, become stable, then, the crystal shape will consist of rough regions 
that have widely separated bunches of size Uf,. The projected free energy for such a surface can be expressed using 
Eq. ( p^ ) as 

/(s)=70+^7i|.s|+%k|^ + O(|sh, (45) 

where the contribution 0(|s|^) comes from the short-ranged attractive interaction between the bunches. In the above 
equation, we have introduced quantities ?/i — [ri(T) + f{nb,T)] /uz and 773 = 7r^(fcBT)^A^/(6n^Ea^). The crystal 
shape z{x) is given by 

Az(a;) = /(77)|„=_A., (46) 
where /(r/) = mins[f{s) — rjs] is the Legendre transform of f{s) to the conjugate field 77, an A is parameter fixing the 



overall size of the crystal |15|. The minimization leads to the result 

z{x) — 0, X < Xo, 
2VX 



X > Xo, (47) 



where xo = "Ui/X. Notice that the crystal the rough regions joins the flat region in a "continuous" with the well known 
"3/2" exponent. The above expression for crystal shapes holds in regions were there is a bunch of a stable size, so 
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that free energy of the surface is obtained by perturbative methods. In the critical region, where we do not have an 
analytical expression of free energy, we will not be able to write down a simple expression like Eq. (^^. However, from 
the general arguments presented in the previous section and the shape of the free energy shown in Fig. , we expect 
the shapes to evolve smoothly as the stable bunch size changes. Finally we note that at very high temperatures, where 
the stable bunch size becomes one, the free energy reduces to the standard free fermion energy and one recovers the 
expression for the crystal shape obtained in Rcf. Q . 



B. 3-Dimensional Crystals 



We now turn our attention from the cylindrical crystals to a 3-dimensional crystal. In this case the crystal has 
flat facets that are joined to the rough regions as shown in Fig. (p^. Once again we will focus on temperatures 
such that the rough regions have stable bunches of size ri},. The orientation-dependent free energy /(p), where 
{dz/dx,dz/dy) — —\p\{cos(j), sincj)), now takes the form 



/(IpI) = 70 + m(0)|p| + mWlpl' + o(IpI^), 



(48) 



where the coefficients rji and 773 are functions of the orientation 0. Note that (p is the angle between the y— axis and 
the average direction of steps at a given point on the surface. In the limit |p| ^ 0, is the angle between the j/— axis 
and the direction of the facet edge. An important quantity that characterizes the three dimensional crystal shape 
is the jump in gaussian curvature at the point where the flat facet joins the rough region. The Gaussian curvature 
K{x,y) of the crystal at any general point x = {x,y) can be written as pq] 



K{x,y) = 



A2(l + |p|2)- 



det [dy{\v\)/dp,dp,] 



p(-Ax) 



from which the Gaussian curvature in at the facet edge (|p| — s- 0) can be cast in the form 

\2 



[6773(</') (771 (</>)+ ^1(0))] 



(49) 



(50) 



Using the fact that 771(1/)) +r]^ {(p) ~ nhJ^{(p)/az and 773(1/)) = 7r^(fcBT)^A^((/))/(67i^I](0)a^), the bunch stiffness, wc find 
that the jump in Gaussian curvature at the facet edge, AK, is satisfies the relation 



kBTVAK 
A 



3/2 



(51) 



We see that the for ansisotropic crystals, where the elastic interactions between the steps depend on their orientation, 
the jump in Gaussian curvature depends on the angle (j). In the limit of high temperatures when A{, — *■ 1 and Uf, — 1, 
we recover the universal relation first derived by Akutsu et al ||26j. In the case when U = 0, so that nt — 1, Eq. ( [5l| ) 
generalizes the expressions that Saam obtained for certain specific values of g, namely, 5 = 0, —1/2 and 2. Note 
that the jump in curvature in addition to being non-universal depends not only on the strength of the long range 
interaction, but also on the strength of the short-range attraction through n;,. We note that in temperature window 
around which the bunch size decreases, we expect the formula derived for the rih bunch to evolve to the corresponding 
expression for the the phase with reduced bunch size. However, in the absence of the exact form of free energy in the 
critical region, we cannot study the details of this evolution. 



VII. COMPARISON OF OUR THEORY WITH CONTINUUM THEORIES 



Recently the effect of competing short range and long range interactions on the finite temperature phase transitions 
on vicinal surfaces was considered by Lassig j|] and Bhattacharjee They studied a continuum version of the 
model considered in this paper, where the fermions interact via a short range attractive part —u of a range a and an 
inverse square repulsive long range interaction of strength g. In the continuum theory, all lengths are measured in 
units of a, which can set to unity without loss of generality. They analyzed the phase transition using renormalization 
group(RG) techniques. For N fermions on ring of size L, such that s = N/L 0, they find that for a given g, 

when the short range attraction is sufficiently attractive, the fermions tend to phase separate. That is, for u > ^^^^^ 
pO|, the fermions on the ring condense into a narrow region. With increasing u, the region's fermionic wave function 
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becomes more and more localized, or in other words the steps try to form a single bunch of macroscopic size N. This 
picture is fundamentally different from our theory. In our picture, the steps do not phase separate to form a single 
macroscopic bunch. Instead, they form bunches of a finite size. The bunches then fill the space and interact with 
each other with renormalized interaction parameters as discussed in the text. With increasing temperature the steps 
undergo a series of "peeling" transitions where the bunch size progressively decreases by one. As mentioned earlier, 
the size of a bunch and peeling transitions are lattice effects which are not present in the continuum model. In the 
continuum model, steps (bosons) are treated as point objects. An infinite number of these points can thus collapse to 
a space of arbitrarily small size in order to maximize the short-range attraction. Therefore, for the physical region we 
are concerned with, namely when the range of the short-range attraction is comparable to the "size" of the bosons, 
the continuum model is not valid. Due to the underlying crystal structure, such a region may well be more relevant 
for describing the experimental situation of surfaces. 

For a fixed strength of the interaction, the continuum theory only captures the onset of the first bunching instability 
as the temperature is decreased. In our model, the first bunching instability that is encountered on cooling from the 
high temperature phase is the point where 2-bunches first appear. For this transition, we find quantitative agreement 
with the continuum RG theory. In particular, the width of the critical region scales with density like ATc ~ 5^1+23 
for 5 < 3/2, which is precisely the relation that we obtain in Eq. (js^). For g > 3/2, the continuum theory fails to 
make any predictions because of a singularity in the fermionic wave functions. However the natural cutoff present in 
the lattice model, prevents such a problem in the lattice model. 

Our theory of step unbinding and the continuum theory also give very different results for the crystal shapes, as 
illustrated in Fig. (p^. Since the steps phase separate in the continuum theory, there is a slope discontinuity at the 
place where the facet of the crystal joins the rough part at low temperatures. This discontinuity in slope decreases 
continuously and vanishes at the tricritical point. Above the tricritical temperature, the rough part joins the facet 
with the 3/2 exponent. In our model the rough part joins the facet with the 3/2 exponent at all temperatures. As 
described earlier, the sizes of the step bunches in the rough regions decrease with increasing temperature. 



VIII. COMPARISON OF THEORY WITH EXPERIMENTS ON SILICON SURFACES 

As noted in the introduction, one of the objectives of the this work is to provide an explanation for the experimentally 
observed exponents and for the physics of bunching phenomenona. In principle, to accomplish this one would have 
to carry out the calculations detailed in the text using interaction potentials between steps derived from electronic 
structure or other quantum mechanical calculations. Since such calculations are not available at this time, we use 
the simple model with an attractive interaction restricted to only near neighbor sites to explain the experimental 
observations. 



A. Experiments of Song et.al 

We first consider the experiments of Song et.al where x-ray scattering studies were conducted on Si(113) 

miscut from the [113] direction towards [001]. It was seen that depending on the slope of miscut surface s, above a 
temperature T'c(s), the surface is uniformly stepped, while the surface consists of bunched steps below Tc(s). The 
number of steps in a bunch was seen to be about 22. From the measurements it was inferred that s ~ |T'c(0) — Tc{s)\^ , 
with P = 0.42 ± .1. Earlier work based on simple mean-field theory M, predicted this exponent to be /? = 1.0. 



Recent renormalization group calculations 
attractive and repulsive interaction strengt 



predict that the exponent depends on the ratio of the strengths of the 
;hs. In the region where these calculations are valid, the exponent P > 0.5. 
The exponent /3 = 0.5 corresponds to a special point for 1-D systems interacting via the inverse square potential, a 
point beyond which the renormalization group methods fail. The continuum theories therefore require a particular 
value for U/G to obtain agreement with experiments. Within our picture we identify the experimentally observed 
transition as a 1-bunch to 22-bunch transition. While this transition has not been explicitly studied in this paper, we 
argue that the features of such a transition can be understood based on results obtained in Sec. IV. First note that 
like the 1-3 transition that we observe for u/g > 2.4, a 1-22 transition can easily be expected in a system with a short 
ranged potential with a more complex structure. Even in a potential with next to near neighbor interaction i.e. rig = 2 
and with U2 > 0.25, at T = 0, one directly goes from a one bunch to a 3-bunch bypassing the 2-bunch. For potentials 
with more complicated features, it is conceivable that one may go directly from a 1-bunch to a phase with many 
steps in a bunch even at T = 0. In these systems, one may easily expect a 1-22 transition, just as we observed a 1-3 
transition in a simple model. Furthermore, the arguments presented in Sec. IV. on the shift in transition temperature 
very well apply to the 1-22 transition provided that the free energy curves /(22,T) and /(1,T) intersect with a finite 
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slope. If this is the case, we find that the shift in transition temperature goes like from which it follows that /3 = 0.5, 
in agreement with the experimentally observed exponent. Another key point is that when the bunch interactions are 
taken into account, we have argued that the transitions that are first order in the EDL become continuous. This 
means that the the exponent (3 describes a curve of continuous 1-22 bunch transitions rather than a curve of first 
order transitions associated with a tricritical point. 

From an analysis of the x-ray scattering amplitudes Song et.al |^ were able to extract the temperature dependence 
of the stiffness of the miscut surfaces. Within our bunching picture the stiffness of a bunched vicinal surface can 
be expressed as |l5| ] 7 = 1 7j_ , where the stiffness along the direction of the steps can be written as 7| | = Sf,/0 
while the stiffness in the direction perpendicular to the steps is expressed in terms of the surface energy of the miscut 
surface 7(6') as 7j^ = 7(0) -I- ^"{9). Using Eq. (^), we can write the scaled surface stiffness, 7 in a stable rif, bunch 
phase as 



Using the value of scaled stiffness measured |^ for the 2.1" miscut surface, 7 = 1.5 at (T — Tc)/Tc ~ 0.15 (here Tc is 
the temperature at which the 22-bunch phase emerges) to represent the stiffness of the ni, ~ 1 phase as the critical 
region is approached and using Eq. (p^, we find g{Tc) « 3/2. Using this value of g and letting n?, — 22, we find the 
scaled stiffness of the bunched phase at {Tc — T)/Tc = .15 to be 1.3. This is closer to the experimentally observed value 
of 2.2 than the value of 0.87 predicted by earlier mean-field studies We note that there are some uncertainties in 
the the experimental results for stiffness in the 22-bunch phase and at low temperatures, which might be the cause of 
the discrepancy between our results and the experimental observations. In Ref. 1^, the authors derive an exponent for 
the temperature dependence of the stiffness by identifying a spinoidal temperature. Within our picture the transition 
is continuous and a spinoidal temperature does not exist. 

Finally, we note that Song et al have studied equilibrating bunches and have sucessfuUy explained the limiting 
bunch size in terms of the Marchenko groove picture discussed in Section 111. A above. 



B. Experiments of Sudoh et. al 

We now turn to the recent experiments of Sudoh et.al [Bj where the morphology of Si(113) surfaces, miscut towards 
a low symmetry azimuth which is 33'^ away from the [110] direction to the [332] direction was studied using scanning 
tunneling microscopy. The authors compared surfaces quenched after annealing at various temperatures in the range 
600-1000°C. Annealing at temperatures above 720*^0 resulted in single, double, and triple-height steps, while annealing 
at a temperatures in the range 690- 720^^0 the average terrace spacing increased due to presence of double, triple, and 
quadruple-layer steps. The size of the terraces saturate during annealing, indicating that the surface had reached 
local thermal equilibrium. STM images at 710*^0 show a predominance of double steps. The authors also find that 
the stiffness of the a step is proportional to its height, with an additional stabilization of double height steps. 

Most of the features of the above observations can be understood using the model considered in this paper, in which 
the attractive interaction is restricted to the nearest neighbor site . We conjecture that the temperature around 720°C 
corresponds approximately to u/g « 2.35 and l/g « 0.5, where the stable phase is the 2-bunch phase. This point is 
marked as "Su" in Fig. (^. Note that both the 1-2 line and 2-3 lines are very close to this point, while the 3-4 line is 
also nearby. This means the free energies of the 1,2 and 3 bunches are very close to each other with the free energy 
of the 2-bunch being the lowest, explaining its stability. Though the 2-bunch is stable in the EDL, the point Su is 
very likely to be in the critical region, where the 1 and 3 bunch phases are stabilized by bunch interactions. Further 
decreasing the temperature results in the emergence of the 4-bunch phase, which can be understood by noting that 
the system then moves closer to the 3-4 line resulting in the stabilization of the 4-bunch due to bunch interactions. 
We also point out our arguments that the stiffness of the bunched steps should be proportional to the number of steps 
in the bunch is also borne out in the experiments. Also, note that this vicinal surface oriented in a direction different 
from that of Song et.al. As a result we expect different strengths of attractive and repulsive step interactions leading 
to different bunch sizes in the two experiments. 

Step bunches with rib in the range 1-4 have also been seen on miscut Si(113) by other workers, specifically van 
Dijken et al ]pl] and Zhu et al [p9|. Equilibrium does not appear to have been achieved in these cases. 
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IX. DISCUSSION AND FUTURE DIRECTIONS 



In this paper we have explored the phase transitions arising in vicinal surfaces with competing long-range and short 
range interactions. We showed that depending on the strength of the interactions and the temperature, the steps on 
the surface rearranged themselves into bunches. Using a tractable model for step interactions, we showed that the 
bunch size increased with decreasing temperature through a series of bunching transitions. At T = 0, we showed 
that the bunch phases predicted by our model correspond to well known periodic groove structure first predicted 
by Marchenko. Our theory can then be considered as an extension of the Marchenko theory to finite temperatures, 
where changes in the size and number of steps in a groove are brought about by the bunching transitions. The 
implications to experiments and the physics of crystal shapes was pointed out. We found that the exponent relating 
the dependence of the transition temperature on the miscut angle, /3, expected from our theory to be consistent with 
the experimentally measured exponent and that the observed features of step stiffness was also in agreement with 
our theory. We now point out some directions for future theoretical and experimental research that can provide more 
insights on equilibrium bunching on vicinal surfaces. 

On the theoretical/computational side, it is very important to have an accurate determination of the step interaction 
potential. This can be done through quantum calculations or atomistic simulations using reliable empirical inter- 
atomic potentials. Once the step interactions are known, the methods developed in this paper can be used to study 
the step bunching transitions. On the experimental side more data on surfaces with competing interactions can shed 
light on the role of attractive interactions on equilibrium bunching properties. In particular for vicinal surfaces in 
a bunched phase, the crystal shape will exhibit a Pokrovsky-Talapov transition as the facet is approached. This is 
also true within the Marchenko picture. An experimental check of this feature would be valuable. It would also be 
useful to determine the exponent (3 for a number of different surfaces that have attractive interactions between steps. 
In contrast to the recent RG calculations that require a particular ratio of the attractive to repulsive interaction 
strengths to obtain (3 = 0.5, our analysis shows that this value of the exponent is robust, applying to all the bunching 
transitions that are first order in the EDL. 

A detailed study of the step phases as a function of orientation of the silicon surfaces away from the [113] direction 
would be most useful in order to elucidate the connection between the experiments of Sudoh et al and those of Song et 
al. Note that miscuts such as those in the former experiments introduce kinks in steps on the Si(113) surface. For small 
kink densities, the model parameters G, J7, and rj should have corrections proportional to the kink density. The phase 
boundaries in Fig. ^ will then shift linearly with kink density, providing a means of tuning the model parameters. 
For cases such as that marked by Su in Fig. (|^), small changes in the kink density would lead to pronounced effects 
due to the close proximity of phase boundaries. 

The quantum n-bunch, or n-mer, phases discovered in this work are highly interesting in their own right, independent 
of their realization in crystal surface physics. While we expect these phases to be Luttinger liquids well away from 
transitions from one bunch phase to another, the picture may change in the crossover region where there is a mixing 
of the long-wavelength modes associated with phonon-like oscillations of the bunches and the modes internal to the 
individual bunches. In effect, new many-body phases may emerge in the critical region. A detailed study of the 
crossover region will shed light on some these new physical phenomenon. 
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APPENDIX A: FREE ENERGY OF AN INCLINED SINGLE STEP IN THE TSK MODEL 

In this appendix we derive the free energy of a single step placed on a rectangular lattice with lattice constants 
{ax,ay)- The step runs on the average at an angle a with respect to the y-axis (The straight step has a — 0.). The 
slope of the step with respect to the y-axis is then I = tana. If the free energy per unit length of the step as a 
function of s is /(s), then in the expansion for small s 

m = fo + l^s' + ■■■, (Al) 

S is the step stiffness. To compute /(s), it is easier to first compute the Legendre-transformed free energy 
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fiv) = fiS) - fis, (A2) 

where ?} is a field coupling to the slope s ||l5[] . We assume that in going from one row to the next, the step can follow 
only three paths. It can go straight to the next row up with energy eo, to the left one column and then up one row 
with total energy eg + F — 77, or to the right one column and then up one row with total energy eg + T + ry. The 
partition function for the step is readily found, yielding the free energy 

ayfifj) = eo - fcsrin[l + e-'^^i'-") + e~^^^+'-'^] (A3) 



Combining Eqs.(Al- A3) yields (noting that e '^'^ ^ 1 for this model to apply), we obtain Eq. (03) of the text. 



APPENDIX B: QUANTUM MONTE CARLO METHODS 

Variational Monte Carlo method: This is the simplest way of performing a quantum Monte Carlo simulation. 
Here, our aim is to have a trial wave function that has a set of parameters that can be varied. One then optimizes 
these parameters so as to obtain the lowest possible energy. By means of a Monte Carlo simulation, one can compute 
the energy of a wave-function as a function of the parameters. 

We start with writing the trial state jV'T > in terms of the configurations of the one dimensional lattice system 

\^pT >=^\R>< RllJjT >, (Bl) 
R 

where the coefficients < i?|'0T > depend on the set of variational parameters Pi. Here, |i? > is a state in configuration 
space. For example, for N bosons on L sites a configuration could be a state with bosons on lattice sites [ii,i2--iN] 
such that 1 < ife < L for 1 < fc < A^. The expectation value of the energy in the state \ipT > can be expressed as 

< VtIHIV't > 

-c^T = ; — n 

< Vt\wt > 

J2r < ^t\R >< R\^T > ' 

where the local energy in configuration R is defined as 



(B2) 



<i^T\n\R> 

^'^^^ - <^Pt\R> ■ ^^^^ 



The sum in Eq. (B2), in general impossible to evaluate analytically, can be evaluated by Monte Carlo. One considers 



the following expression: 



where one interprets 



i?T-^£;L(i?)p(i?), (B4) 

R 



(R\ <^T\R><R\ijT> 

'^""^ ^ Y.r<^t\R><RHt> ^^'^ 

as a probability density function. If one then samples M configurations distributed according to p{R) then the exact 
expression for Et can be approximated by 

p _ Y.'rEl{R) 

EvMc — J y^») 

where the prime is used to represent sum over the configurations. In the limit of large M, Evmc will converge to 
Et- The average Evmc, of course, is subject to statistical fluctuations, and one can easily evaluate the statistical 
error by SE = cje/VM. The variance of the local energy, as, is defined as 



= \/<E^>- <E >2, (B7) 
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where < E >= Evmc and < E"^ >— ' — • It is known pi| that as, which has a lower bound of 0, is a better 
quantity to optimize than Evmc- We use both quantities as indicators in our search of variational parameters. 

For each set of variational parameters, we use the Metropolis algorithm to generate a set of M configurations 
distributed according to p{R)- We start with a configuration of M walkers each of which has N bosons randomly 
distributed on a lattice of one dimensional ring of size L. Each walker then undergoes a "stochastic walk" according 
to the following procedure: (1) for each boson we pick a near neighbor site randomly and (2) if the site is occupied 
the boson is not moved, but if the site is empty the boson is moved with probability q given by 



min 



V't \ ne'w 



(B8) 



where ipTlnewioid) is (V'tI^) (RIiPt)] evaluated at Rnew and Roid- This process is repeated sufficient number of times 
so that one ends up with a set of walkers distributed according to p{R). 

Green's Function Monte Carlo method: This is in principle an exact way to calculate the ground state properties 
of Bose systems. Its required computer time scales algebraically (as opposed to exponentially in exact diagonalization) 
with system size, thus allowing exact calculations for large systems. Here, one combines random sampling with a 
simple method to project out the ground state from an arbitrary initial state. 
We start with the operator 

T = C-n, (B9) 

where the constant C, whose choice we discuss below, is to ensure that all matrix elements {R'\J-'\R) are non- negative. 
The basic premise of GFMC is that repeated application of on an essentially arbitrary initial state will project out 
the ground state. That is, if an initial state > has any overlap with the ground state lipo) of Ti, the process 

|^(«) >^^"|^(o) > (BIO) 
will lead to IV'o) at large n. GFMC is a method to realize the above process by a Monte Carlo random walk. 



In order to improve efficiency of the random walk process, one more mathematical manipulation of Eq. ( BIO ) is 
necessary. This is done by introducing an operator whose matrix elements are related to those of by a similarity 
transformation: 

f{R\R) = (V^Tli?') {R'\T\R)/{-iPt\R). (Bll) 



The stochastic realization of Eq. (BIO) is actually with JF instead of J-. While mathematically equivalent, the use of 
can significantly reduce the fluctuation of the Monte Carlo process if \iPt) is a reasonable approximation of |V'o)- 
This is known as importance sampling [ p3[ . 

The program is then to start with a set of walkers distributed according to < iPt\R >< Rl'fpT >■ These walkers 
undergo random walks in i?-space. For each walker, denoted by R, taking a step in the random walk means randomly 
selecting and moving to a new position R' with probability T{R\ R)/ '^ji> ^{R', R)- Note that, due to the structure 
of Ti, there are only 0{N) possible i?"s and the corresponding probabilities can be easily computed. Because the 
overall normalization J-"(i?', i?) is not a constant, each walker carries an overall weight which fluctuates, or a 
branching scheme is introduced to allow the total number of walkers to fluctuate, or both. 
The ground state energy is given exactly by the so-called mixed estimate 

< iJT\n\tPo > 
-^0 ^ ^ — r, 

^ T.rEl{R)<^Pt\R><RHo > , . 

where E]^{R) is defined in Eq. (p3|). After a sufficient number of steps, the walkers are distributed according to 
< VtI^ >< ^IV'o >• can therefore be computed from such walkers as the (weighted) average of El{R) with 
respect to walker positions R, similar to Eq. (|B^). We denote this Monte Carlo estimate of Eq by Eqpmc- 

Expectation values of operators other than Ti can also be computed from the mixed estimate. However, if the 
operator does not commute with ?i, this estimate is not exact. As we mention in Section this is an important 
distinction which requires careful analysis of the results. There exist ways to improve upon the mixed estimate and 
to possibly extract exact estimates of expectation values pl,p4], but we will not discuss them here. 
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Now we describe our choice of C. Since we are dealing with a finite number of hard core bosons on a lattice of 
finite size, the eigenstates of Ti. are bounded both from above and below. If Emax is the highest eigenvalue of 7i, we 
choose C such that C > Emax- This is easily achieved for the problem at hand. On the lattice we put all bosons next 
to one another and compute the potential energy arising from the long range interactions alone. For N bosons, we 
choose C as 

C = 2N + gY,^^^. (B13) 

i=l 

This clearly satisfies C > Emax, since hopping and the short range attraction can only lower this energy. 
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FIGURE CAPTIONS 



FIG. 1. Geometry for the TSK model. The steps between terraces are indicated by bonds on the lattice. They run parallel 
to the y-axis on an average. 

FIG. 2. T=0 phase diagram of the vicinal surface plotted as a function of the ratio of the near neighbor attractive interaction 
U and repulsive inverse square interaction G. The dots separate regions with stable bunch sizes that differ by one (The bunch 
sizes are given indicated in the figure). 

FIG. 3. Plot of the difference in energy per boson in a nt bunch and the energy of a 1-bunch. At the points indicated by 
dots, the bunch size corresponding to the minimum of the energy difference increases by one. 

FIG. 4. Periodic Groove structure proposed by Marchenko. The figure shows the repeating groove structure (length Lg), the 
flat facet (length Li) and the stepped facet (length L2). The case n;, = 3 is shown here for graphical simplicity. Marchenko's 
calculation applies only for nj, large. 

FIG. 5. Energy/step in bunches of size 2(circle), 3(square) and 4(triangle) for u/g = 1.65, plotted as a function of g. The 
energies are obtained by exact diagolization of the quantum bunch problem on 12 sites in the case of 3 and 4 bunches and 200 
sites for the 2 bunch. 

FIG. 6. Phase diagram in the extreme dilute limit plotted in the u/g-l/g space. The lines separate regions with stable 
bunch sizes that differ by one. The stable bunch size in each region is indicated in the figure. The 1-2 line (dotted line) lies 
above the 1-3 line (short-dashed line) for u/g < 2.4, while lies beneath it for u/g < 2.4. For u/g > 2.4 the 1-3 line lies in 
between the 1-2 and 2-3 (solid line) lines. This implies that for u/g > 2.4, the 1-bunch phase directly undergoes a transition 
to the 3-bunch phase bypassing the 2-bunch phase. The long-dashed line is the 3-4 line. 

FIG. 7. Comparison of the exact diagonalization calculations and GFMC simulations. The circles indicate the energy per 
step in a 3-bunch confined to 12-sites and for a 2-bunch on 200 sites as found from exact diagonalization, while the triangles 
and squares represent the results for the 2-bunches and 3-bunches on 30 sites and 100 sites, respectively, as found using GFMC. 
The error bars on the GFMC energies are shown in the figure. Lines are drawn through the points to guide the eye. 



FIG. 8. Plot of the 1-3 line and 3-4 line for u/g < 8. 

FIG. 9. Schematic plot of the energy per step (solid curves) when the bunch-bunch interactions are taken into consideration. 
The dotted lines indicate the energy per step in the EDL. In Region 1, the nt-bunch has lower energy while in Region 3 the 
Ub ~ 1 bunch is lower in energy. In the critical region, the energy per step smoothly passes from f{nt,T) to f{nt — 1,T). The 
width of the critical region indicated in the figure is estimated in the text. 

FIG. 10. Schematic depiction of the two body correlation f{r) used in the trial wave- function. The function f{r) takes on a 
value fiair — l has a minimum value of d at r = ro and behaves like 1 — const /r" for large r. 

FIG. 11. The boson pair correlation function g{x) for 10 bosons on 100 sites with u/g — 1.25 for values of g indicated in the 
figure. The points are also marked in Fig. (^^ for comparison. By counting the number of peaks in g{x) and making use of 
the periodic boundary condition the number of bosons in a bunch can be established. For example, in the case g — 20, there 
are 2 peaks at a; = 20 and a; = 40 in addition to a large value of g{x) on the near neighbor site. Invoking the periodic boundary 
condition, one sees that if there are 2 bosons in a bunch, the four peaks at a:: = ±20, x = ±40 account for 8 bosons, which 
along with the boson on the near neighbor site and the boson at the origin add up to 10 bosons. For g — 20 and g = 15 we 
have a well defined 2-bunch phase, while there is a well defined 1-bunch phase for g — 10. The figure shows that g = 12 is in 
the critical region. We illustrate the statistical error bars only for the case g = 20, since in all the other cases they have similar 
magnitudes. 
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FIG. 12. Boson pair correlation function (Eq. (|l3|)) for 12 bosons on 120 sites with u/g = 1.65 for values of g indicated in 
the figure. The points are also marked in Fig. (|l^) for comparison. Counting the number of peaks as we did in Fig. (^l|) gives 
the number of bosons in a bunch. To account for the 3-bunches, note that when g = 15 there are 2 peaks at a; = 30 and 
a; = 60 in addition to large values of g{x) on the near-neighbor and next to near-neighbor sites. Invoking the periodic boundary 
condition, we see that the peaks at a:: = ±30 account for 6 bosons while the peak at a:: = 60 accounts for 3 bosons which along 
with the boson at the origin and the ones on two adjacent sites account for all 12 bosons. For g = 15 and g = 7.5 we have a 
well defined 2-bunch phase, while there is a well defined 2-bunch phase for g = 5 and 1-bunch phases for p = 4 and 3 = 3. The 
figure shows that g = 6 is in the critical region. The statistical error bars for all values of g have about the same magnitude as 
the case g = 15 where they have been explicitly displayed. 

FIG. 13. Boson pair correlation function (Eq. (|43| )) for 12 bosons on 120 sites with u/g = 2.0 for values of g indicated in 
the figure. The points are also marked in Fig. (hTfor comparison. By counting the number of peaks in g{x) as we did in 
the previous figures, we see that for g = 7.0, g = 3.5, g = 2.85 and g — 2.25 we have well defined 4,3,2 and 1-bunch phases 
respectively, while the points g — i and g — 2.96 are in the critical region. The error bars have only been shown for the case 
g = 7, since in all the other cases they have similar magnitudes. 

FIG. 14. Comparison of the phase diagrams obtained in the EDL using exact diagonalization and GFMC simulations. For 
u/g = 1.25,1.65,2.0 the pair correlation functions for the marked values of g are given in Fig. (|lll)-(|l^). For each value 
of g we also indicate the bunch size (lb, 2b, 3b and 4b correspond to 1,2,3 and 4 bunches respectively, while cr refers to the 
critical region) obtained from GFMC simulations. Note that the GFMC results and the exact diagonalization results do not 
agree exactly since in the dilute limit the bunch interactions can alter the nature of the phases especially close to the phase 
boundaries. 

FIG. 15. Top view of a 3-d crystal. The flat facets and the rough regions are shown in the figure. The steps in the rough 
region run parallel to the direction inclined to the vertical by an angle (j) as shown, for a point at the facet edge, in the flgure. 
The z-axis points outside the plane of the paper as indicated in the flgure. 

FIG. 16. Comparison of the temperature dependence of crystal shapes in lattice and continuum theories. In the lattice 
theory the rough region always joins the facet continuously with the standard "3/2" exponent. In this theory, the rough 
regions consists of steps in bunches, whose sizes decrease with increasing temperature. Within the continuum theory, the rough 
regions joins the facet with a discontinuous slope. This discontinuity decreases with increasing temperature and vanishes at 
the tricritical point. 
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TABLES 



TABLE I. Comparison of the exact diagonalization and GFMC results for 3 bosons confined on 12 sites with u/g = 1.65. 
Eexad denotes the energy per boson obtained by exact diagonalization, Evar is the variational energy per boson on 12 sites and 
E(3FMC represents the GFMC result. The number bracket, corresponding to the GFMC energies, indicates the error on the 
last decimal place. In all the cases a variational wavefunction with ro = 3, /i = 1.14, d = 0.5 and a = 0.6 was chosen. We also 
give the GFMC energies ,E^pi^jQ andEQ^pjijQ , for 3 bosons occupying 30 and 100 sites respectively. 

_g Eexad EvMC E}?pf^.[Q 

6^5 -0.44902 -0.34(3) -0.4489(2) -0.4490(4) -0.4492(7) 

6.9 -0.57562 -0.47(3) -0.5756(1) -0.5756(4) -0.5753(6) 

7.3 -0.70425 -0.63(4) -0.7042(1) -0.7043(4) -0.7046(6) 



^GFMC 



TABLE IL Energies for u/g = 1.25 and values of g shown in the table, computed for 10 bosons on 100 sites. The variational 
energy is denoted by Evar, the GFMC energy by Ecfmc and Eapprox is the approximate energy computed using Eq. (p2[). 
The variational wave function parametrized by ro, /i, d and a (see Eq. (El) is also tabulated. 



g 


ro 


,/i 


d 


Q 


EvMC 


Eqfmc 


Eapprox 


10.0 


2.0 


1.05 


0.58 


1.5 


3.63(4) 


2.312(1) 


3.071 


12.0 


2.0 


1.06 


0.58 


1.5 


3.76(5) 


2.221(1) 


3.553 


15.0 


2.0 


1.06 


0.30 


0.5 


2.37(6) 


-0.479(1) 


-0.424 


20.0 


2.0 


1.14 


0.18 


0.5 


-0.10(1) 


-6.701(2) 


-6.223 



TABLE IIL Energies for u/g = 1.65 and values of g shown in the table, computed for 12 bosons on 120 sites. The symbols 
have same meanings as in Table. II 



g 


ro 


/i 


d 


a 


EvMC 


Egfmc 


Eapprox 


3.0 


4.0 


1.02 


0.7 


1.0 


1.99(2) 


1.139(1) 


1.311 


4.0 


4.0 


1.02 


0.3 


0.4 


1.85(4) 


1.032(1) 


1.579 


5.0 


4.0 


1.02 


0.4 


0.4 


0.40(3) 


-0.620(1) 


-0.486 


6.0 


4.0 


1.14 


0.1 


0.1 


-1.68(9) 


-3.340(3) 


-3.363 


7.5 


4.0 


1.14 


0.1 


0.1 


-6.83(7) 


-8.870(7) 


-8.742 


15.0 


4.0 


1.16 


0.1 


0.1 


-33.77(3) 


-39.075(3) 


-38.814 



TABLE IV. Energies for u/g = 2.0 and values of g shown in the table, computed for 12 bosons on 120 sites. The symbols 
have same meanings as in Table. II 



g 


ro 


/i 


d 


Q 


EvMC 


Egfmc 


Eapprox 


2.25 


3.0 


1.04 


0.7 


0.7 


1.11(1) 


0.807(1) 


1.104 


2.85 


3.0 


1.06 


0.5 


0.1 


0.89(1) 


0.051(5) 


0.119 


2.96 


4.0 


1.02 


0.1 


0.1 


0.34(2) 


-0.301(1) 


-0.209 


3.5 


4.0 


1.02 


0.1 


0.1 


-2.31(3) 


-3.181(1) 


-3.160 


4.0 


4.0 


1.06 


0.1 


0.1 


-4.61(1) 


-6.165(5) 


-6.144 


7.0 


4.0 


1.06 


0.1 


0.1 


-24.21(3) 


-26.779(1) 


-26.665 
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